GLMM Maternal Mortality Federally- Summer 2025

This is a Report Template

Slides: slides.html ( Go to slides.qmd to edit)

Important

Remember: Your goal is to make your audience understand and care about your findings. By crafting a compelling story, you can effectively communicate the value of your data science project.

Carefully read this template since it has instructions and tips to writing!

Introduction

Generalized Linear Mixed Models (GLMMs) are a flexible class of statistical models that combine the features of two powerful tools: Generalized Linear Models (GLMs) and Mixed-Effects Models(Agresti 2015). Like GLMs, GLMMs can model non-normal outcome variables, such as binary, count, or proportion data. However, they go a step further by incorporating random effects, which account for variation due to grouping or clustering in the data, correlated observations, and overdispersion. The GLMM method can work with complex data structures such as temporal correlation, spatial correlation , nested data, heterogeneity of variance, and repeated measurement(Zuur et al. 2009).

In practical terms, GLMMs are especially useful when data points are not independent, such as when students are nested within schools, patients are treated within hospitals, or repeated measures are taken from the same subject over time. For example, Thall wrote that issues with longitudinal clinical trial basic count data from repeated measures taken from the same subject over time will have problems detecting comparable between subject outcomes because it can be difficult to determine if outcomes are time dependent or due to treatment groups, thus a general linear mixed model method may be utilized to represent dependence upon each patient, incorporate covariate data, create time as a function, account for variability between patients,and be flexible and tractable (Thall 1988). The random effects help model the correlation within clusters and allow for unobserved heterogeneity—differences that are not captured by the measured covariates.

GLMMs are good for:

• Handling hierarchical or grouped data (e.g., students within classrooms, patients within clinics)(Lee and Nelder 1996)

• Modeling non-normal outcomes, such as:

  • Binary outcomes (using logistic GLMMs)(Wang et al. 2017)

  • Count data (using Poisson or negative binomial GLMMs)(Candy 2000)

  • Proportions or rates (Salinas Ruı́z et al. 2023)

• Improving inference by accounting for both fixed effects (predictors of interest) and random effects (random variation across groups)

• Reducing bias and inflated Type I error rates that can result from ignoring data structure(Thompson et al. 2022)

GLMMs are ideal when your data is both complex in structure and involves non-Gaussian response variables, making them indispensable in fields like medicine, ecology, education, and social sciences. Tawiah et al describes zero-inflated Poisson GLMMs, an extension of Poisson GLMM that allows for overdispersion due to a prevalence of zeros in the data, common in health sector data(Tawiah, Iddi, and Lotsi 2020). The paper compares a Poisson GLM, a zero-inflated Poisson GLM, a Poisson GLMM, and a zero-inflated Poisson GLMM, applied to clustered maternal mortality data. Another paper by Owili et al utilizes a GLMM to investigate the impact of particulate matter on maternal and infant mortality globally (Owili et al. 2020). They use a Poisson link function and take year and country as random effects to account for differences in global data.

We wish to analyze the data of federal maternal mortality deaths via VSRR Provisional Maternal Death Counts and Rates dataset by utilizing a General linear mixed model with Poisson link as it is count data. We wish to see if ethnicity (a fixed effect) has any influence upon maternal death per live births(number of live births per 12 month period or year) count by year(random effect). Like other public health or clinical data there will be issues such as correlated observations and overdispersion but GLMM will be utilized to parse through the noise and determine if indeed there are some patterns of maternal mortality among mothers of differing ethnic ties.

Methods

Math Background

GLMMs can be considered an extension of GLMs, wherein a GLM includes the addition of random effects, or an extension of Linear Mixed Models (LMMs), where a linear model with fixed and random effects is extended for non-normal distributions(Salinas Ruı́z et al. 2023). Let

• \(\mathbf{y}\) be a \(Nx1\) column vector outcome variable

• \(\mathbf{X}\) be a \(Nxp\) matrix for the \(p\) predictor variables

• \(\boldsymbol{\beta}\) be a \(px1\) column vector of the fixed effects coefficients

• \(\mathbf{Z}\) is a \(Nxq\) matrix of the \(q\) random effects

• \(\mathbf{u}\) is a \(qx1\) vector of random effects, and

• \(\boldsymbol{\epsilon}\) is a \(Nx1\) column vector of the residuals

Then the general equation for the model is given by:

\[\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\mathbf{Z}{u}+\boldsymbol{\epsilon}\] (Salinas Ruı́z et al. 2023). The GLMM Model process is that the analyis of variance model or the equation is a vector of linear predictors with of unknown parameters estimates. Each distribution has is its own probability function which we will utilize the Negative Binomial as GLMMs typically include a link function that relates the response variable \(\mathbf{y}\) to a linear predictor, \(\eta\), which excludes the residuals. So then \[\boldsymbol{\eta}=\mathbf{X}\boldsymbol{\beta}+\mathbf{Z}\boldsymbol{\lambda}\] The link function is \(g(\cdot)\), where \[g(E(\mathbf{y}))=\boldsymbol{\eta}\] where \(E(\mathbf{y})\) is the expectation of . The choice of link function depends on the outcome distribution. For this paper our data demonstrates a Negative Binomial distribution for overdispered count data, so we will use a log link function.
\[g(\cdot)=log_e(\cdot)\] In the GLMM model the parameter estimates is solved by reducing the negative log likelihood functions(Salinas Ruı́z et al. 2023). The means or the least square means are derivative of the parameter estimates and are found on the model scale. The link function, negative binomial log link, will convert the mean estimates at the model scale to the data scale. Negative Binomial Distribution: \[ f(y;k,{\mu})=\frac{\Gamma(y+k)}{\Gamma(k)*(y+1)}\left(\frac{k}{\mu+k}\right)^{k}\left(1-\frac{k}{\mu+k}\right)^{y} \] Zuur writes that the negative binomial Distribution has two parameters \({\mu}\) and \(k\) (Zuur et al. 2009). The symbol \({\Gamma}\) is defined as \({\Gamma(y+1)=(y+1)!}\) The Mean of Negative Binomial is given: \(E(Y)= {\mu}\) The Variance of Negative binomial is given; \(Var(Y)= {\mu}+ \left(\frac{\mu^2}{k}\right)\), where second term determines the overdispersion, \(k\) is called the dispersion parameter and indirectly determines overdispersion. If \(k\) is significantly large relative to \({\mu^2}\) then the second term will approximate to zero and a Poisson distribution may as well be used. However, the smaller the \(k\) value the larger the overdispersion may form and then negative binomial is the correct log link to utilize.

GLMM model process is

Assumptions

Before deciding to use a GLMM for our data, we had to check some assumptions (specific to our negative binomial distributed data).

• The response variable and the predictors have a linear relationship within the levels of random effects.
• The response variable is assumed to follow a negative binomial distribution, with \(\sigma^2>\mu\).
• The residuals and random effects are independent.
• The random effects are assumed to be normally distributed, with mean 0 and variance \(\sigma\).

Model Choice

The Negative Binomial distribution’s quadratic mean-variance relationship is ideal for overdispersed count data that the Poisson distribution cannot accommodate. Overdispersion is a given for populations due to their heterogeneity that aggregation processes such as community clusters form with shared traits such as ethnicity. Our data includes community cluster groups such as ethnicity and age groups. The negative binomial model is desirable also due to its easily interpretable dispersion parameter as a measure of aggregation and its tractability due to closed form expression of probability mass function which helps in convenient model inference and estimation. Since our data is longitudinal, a Negative Binomial GLMM is ideal since measurements are taken over time. We incorporate Year as a random effect because it could explain for variation in our model which could not be explained by our fixed effects such as ethnicity or age group.

We perform the analysis with R(R Core Team 2025) and the packages lme4{Bates et al. (2015)} with function glmer to investigate poisson model or glmer.nb() for negative binomial model and the package glmmTMB{McGillycuddy et al. (2025)} with the function, nbinom2(), for negative binomial model. The approximate parameter estimation method uses for glmmTMB uses AGQ or Laplace approximation with Wald Z significance test. The default approximate parameter estimation method for lme4 is Laplace approximation when no quadrature point is used but is pseudo-likelihood method in linearization,and integral approximation and another estimation method is Gauss-Hermite quadrature to approximate the log-likelihood using numerical integration(Salinas Ruı́z et al. 2023).

Densities are seen as counts per volume and can be modeled with NB and offset variable(Zuur et al. 2009). We are investigating maternal deaths per live births and while we can use the variable Maternal Mortality rate, it would be less noisy if we used the offset variable log(Live Births). We can also model without the offset and utilize Akaike information Criterion,AIC, to see which is the better model as the smaller the AIC would indicate a better model. A random effect could be month,year, or month/year. Fixed effects could be age groups, ethnicity, and perhaps year. Our choice of random or fixed effects are dependent on the question we which to investigate.

Analysis and Results

Data Exploration and Visualization

The data we used comes from the National Vital Statistics System, an organization that provides national data about births, deaths, marriages, divorces, and fetal deaths. It is a collaboration between the National Center for Health Statistics (NCHS) (a part of the CDC) and state vital records offices. {create citation https://www.cdc.gov/nchs/nvss/about_nvss.htm} This dataset is the Vital Statistics Rapid Release (VSRR) Provisional Maternal Death Counts and Rates, in the form of a .csv file. The dataset contains monthly death counts and death rates from 2019 to 2024 by race and ethnicity, age and overall. Rates are maternal deaths per 100,000 live births. For this dataset, a maternal death is defined as “the death of a woman while pregnant or within 42 days of termination of pregnancy irrespective of the duration and the site of the pregnancy, from any cause related to or aggravated by the pregnancy or its management, but not from accidental or incidental causes.” {create citation https://www.cdc.gov/nchs/nvss/vsrr/provisional-maternal-deaths-rates.htm}

• Present initial findings and insights through visualizations.

• Highlight unexpected patterns or anomalies.

# loading packages 
library(tidyverse)
library(knitr)
library(ggthemes)
library(ggrepel)
library(dslabs)
library(ggplot2)
library(patchwork) # For combining plots
library(ggpubr)
library(DataExplorer) #For EDA
options(repos = c(CRAN = "https://cloud.r-project.org"))#specify cran mirror
if (!require("lme4")) install.packages("lme4")#if required to install
library(lme4)
##to check if model is overdispersed'
if (!require("performance")) install.packages("performance")#if required to install
library(performance)
##Overdispersion is detected with Poisson model thus we must use negative binomial model
if (!require("glmmTMB")) install.packages("glmmTMB")#if required to install
library(glmmTMB)

df <- read.csv("data/VSRR_Provisional_Maternal_Death_Counts_and_Rates_20250726.csv")#loading data
#variables- Maternal.Deaths--outcome variable
#Live.Births---offset variable since it is maternal deaths per live births
#Subgroup- ethnicity and year of mother must be edited --- these are fixed
#Year.of.Death---random effect variable
#Summary
#renaming columns
df<-df%>% 
  rename(Year= Year.of.Death, Month= Month.of.Death, Maternal_Deaths= Maternal.Deaths, Live_Births= Live.Births, Maternal_Mortality_Rate= Maternal.Mortality.Rate, Time_Period= Time.Period, Month_Ending_Date= Month.Ending.Date, Data_As_Of= Data.As.Of) 
str(df) 

'data.frame':   720 obs. of  12 variables:
 $ Data_As_Of             : chr  "04/06/2025" "04/06/2025" "04/06/2025" "04/06/2025" ...
 $ Jurisdiction           : chr  "United States" "United States" "United States" "United States" ...
 $ Group                  : chr  "Total" "Total" "Total" "Total" ...
 $ Subgroup               : chr  "Total" "Total" "Total" "Total" ...
 $ Year                   : int  2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
 $ Month                  : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Time_Period            : chr  "12 month-ending" "12 month-ending" "12 month-ending" "12 month-ending" ...
 $ Month_Ending_Date      : chr  "01/31/2019" "02/28/2019" "03/31/2019" "04/30/2019" ...
 $ Maternal_Deaths        : int  660 653 657 668 706 720 721 723 729 739 ...
 $ Live_Births            : int  3787776 3783489 3771682 3772235 3767999 3757275 3761912 3758847 3761856 3760121 ...
 $ Maternal_Mortality_Rate: num  17.4 17.3 17.4 17.7 18.7 19.2 19.2 19.2 19.4 19.7 ...
 $ Footnote               : chr  "" "" "" "" ...

filtered_df<- df%>% filter(Maternal_Deaths > 0, Live_Births > 100) #Filtering out rare combinations

test_df <- filtered_df%>%
  select(-Maternal_Mortality_Rate)
  na.omit

function (object, ...) 
UseMethod("na.omit")
<bytecode: 0x2caab68>
<environment: namespace:stats>

#Summary Statistics

deaths_df <-test_df %>%
  mutate(
    Ethnicity = case_when(
      Group == "Race and Hispanic origin" ~ Subgroup,
      TRUE ~ NA_character_
    ),
    Age_Group = case_when(
      Group == "Age" ~ Subgroup,
      TRUE ~ NA_character_
    ),
    Is_Total = Subgroup == "Total"
  ,
  Year = as.factor(Year),
  Month= as.factor(Month),
  Ethnicity=as.factor(Ethnicity),
  Age_Group=as.factor(Age_Group))#mutating groups

deaths_df_ethnic<-deaths_df %>%
  filter(!is.na(Ethnicity))%>%
  mutate(Ethnicity = Subgroup)###filter out is total if ethnicity or age group are na

deaths_df_age <-deaths_df %>%
  filter(!is.na(Age_Group)) %>%
  mutate(Age_Group = Subgroup)


merged_deaths_df <-  bind_rows(deaths_df_ethnic, deaths_df_age)  #merge both df

#setting factor levels
deaths_df3 <- merged_deaths_df %>%
  mutate(
    Ethnicity = factor(Ethnicity, levels = c(
      "Asian, Non-Hispanic", "Black, Non-Hispanic", "White, Non-Hispanic",
      "Hispanic", "American Indian or Alaska Native, Non-Hispanic",
      "Native Hawaiian or Other Pacific Islander, Non-Hispanic", "Unknown"
    )),
    Age_Group = factor(Age_Group, levels = c("Under 25 years", "25-39 years", "40 years and over", "Unknown"))
  )

#mutating any NAs with unknown
deaths_df3<- deaths_df3 %>%
   mutate(
    Ethnicity = replace_na(Ethnicity, "Unknown"),
    Age_Group = replace_na(Age_Group, "Unknown")
  )

# make variable in which we combine Year and Month into a for of a date to create later variable, or use month ending date?
deaths_df3$Date <- as.Date(paste(deaths_df3$Year, deaths_df3$Month, "01", sep = "-"))

# Dobbs was decided June 24, 2022
cutoff <- as.Date("2022-06-24")##when Roe v.Wade was overturned

deaths_df3$Dobbs_Era <- ifelse(deaths_df3$Date < cutoff, "Pre-Dobbs", "Post-Dobbs")
deaths_df3$Dobbs_Era <- factor(deaths_df3$Dobbs_Era, levels = c("Pre-Dobbs", "Post-Dobbs"))


df_age_year <- deaths_df %>%
  filter(!is.na(Age_Group), !is.na(Year)) %>%
  group_by(Year, Age_Group) %>%
  summarise(Maternal_Deaths = sum(Maternal_Deaths, na.rm = TRUE))
df_age_year

# A tibble: 18 × 3
# Groups:   Year [6]
   Year  Age_Group         Maternal_Deaths
   <fct> <fct>                       <int>
 1 2019  25-39 years                  6017
 2 2019  40 years and over            1139
 3 2019  Under 25 years               1320
 4 2020  25-39 years                  6824
 5 2020  40 years and over            1361
 6 2020  Under 25 years               1312
 7 2021  25-39 years                  8338
 8 2021  40 years and over            1975
 9 2021  Under 25 years               1558
10 2022  25-39 years                  9191
11 2022  40 years and over            2035
12 2022  Under 25 years               1796
13 2023  25-39 years                  6262
14 2023  40 years and over            1158
15 2023  Under 25 years               1134
16 2024  25-39 years                  5610
17 2024  40 years and over            1252
18 2024  Under 25 years               1386

df_total_year <-deaths_df%>%
  filter(Group=="Total") %>%
  group_by(Year, Is_Total) %>%
  summarise(Maternal_Deaths = sum(Maternal_Deaths, na.rm = TRUE))
df_total_year

# A tibble: 6 × 3
# Groups:   Year [6]
  Year  Is_Total Maternal_Deaths
  <fct> <lgl>              <int>
1 2019  TRUE                8482
2 2020  TRUE                9503
3 2021  TRUE               11871
4 2022  TRUE               13022
5 2023  TRUE                8554
6 2024  TRUE                8248

#deaths_cleaned <- deaths_df3 %>%
#  filter(Age_Group != "Unknown",Ethnicity != "Unknown") %>%
 # droplevels()

#deaths_df2<-deaths_df %>%
#  filter(!(Is_Total & (is.na(Ethnicity) | is.na(Age_Group))))###filter out is total if ethnicity or age when total is true

The below chart gives an idea of the missing values in the dataset. It seems like there are many NAs in ethnicity and age group, certainly more than is acceptable. The dataset indicates that Maternal Deaths between 1 and 9 are suppressed for privacy reasons. With more investigation, it appears that the racial/ethnic group “Native Hawaiian or Other Pacific Islander, Non-Hispanic” has NA for all but 2 Maternal Deaths values. We decided to omit this subgroup entirely. There are missing Maternal Mortality Rate entries in “American Indian or Alaska Native, Non-Hispanic”, which will not be an issue as we have decided not to use rate in our model ( though we will examine it during our Exploratory Data Analysis). Finally, there are six NA Maternal Deaths also in the “American Indian or Alaska Native, Non-Hispanic” category, a number that is reasonable to omit.

plot_missing(deaths_df)

plot_missing(test_df)

Since our data is longitudinal, it is easy to visualize with graphs. Below is a bar chart that shows the distibution of Maternal Deaths by racial/ethnic group. The counts indicated the sum of the Maternal Deaths for each group over the span of the data (2019-2024).

#Data Exploration Graph 1
colors=c('yellow','orange', 'blue', 'red', 'green')
df_ethnicity <- deaths_df %>%
    filter(!is.na(Ethnicity)) %>%
    group_by(Ethnicity) %>%
    summarise(Maternal_Deaths = sum(Maternal_Deaths, na.rm = TRUE))
barplot(df_ethnicity$Maternal_Deaths,names.arg=c('Am In/Ala', 'Asian, NH','Black, NH','Hispanic','White,NH'),main="Maternal Mortality Totals 2019-2024 by Ethnicity",col=colors,xlab='Ethnicity',ylab='Maternal Deaths')

This graph shows the total Maternal Deaths for each year. The counts for year are the sum of the monthly data.

#Data Exploration Graph 2
df_tot_year_ts<-ts(df_total_year$Maternal_Deaths, start = c(2019), frequency = 1)
plot.ts(df_tot_year_ts,main='Yearly Total Maternal Deaths 2019-2024',xlab='Year',ylab='Maternal Deaths',col='navy')

df_tot_year_rate <- filtered_df%>%
  filter(Group=='Total')%>%
  select(Maternal_Mortality_Rate)
df_tot_year_rate_ts<-ts(df_tot_year_rate$Maternal_Mortality_Rate, start = c(2019,1), frequency = 12)
plot.ts(df_tot_year_rate_ts,main='Monthly Total Maternal Mortality Rate 2019-2024',xlab='Year',ylab='Maternal Deaths per 100,000 Live Births',col='navy')

total_df <- deaths_df%>%
  filter(Group=='Total')
total_ts <- ts(total_df$Maternal_Deaths,start=c(2019,1),frequency=12)
plot.ts(total_ts,main='Monthly Total Maternal Deaths 2019-2024',ylab='Maternal Deaths',xlab='Time',col='navy')  

His <- filtered_df%>%
  filter(Subgroup=='Hispanic')%>%
  select(Maternal_Deaths)%>%
  rename(Hispanic=Maternal_Deaths)
bl <- filtered_df%>%
  filter(Subgroup=='Black, Non-Hispanic')%>%
  select(Maternal_Deaths)%>%
  rename('Black_Non-Hispanic' = Maternal_Deaths)
wh <- filtered_df%>%
  filter(Subgroup=='White, Non-Hispanic')%>%
  select(Maternal_Deaths)%>%
  rename('White_Non-Hispanic'=Maternal_Deaths)
as <- filtered_df%>%
  filter(Subgroup=='Asian, Non-Hispanic')%>%
  select(Maternal_Deaths)%>%
  rename('Asian_Non-Hispanic'=Maternal_Deaths)
aina <- df%>%
  filter(Subgroup=='American Indian or Alaska Native, Non-Hispanic')%>%
  select(Maternal_Deaths)%>% 
  rename('American_Indian_or_Alaska_Native, Non-Hispanic'=Maternal_Deaths)
ethn_ts <- bind_cols(His,bl,wh,as,aina) 
ethn_ts <- ts(ethn_ts, start=c(2019,1), frequency = 12)
colors=c('red','blue','green','orange','yellow')
plot.ts(ethn_ts, plot.type='single',main='Monthly Maternal Deaths by Ethnicity 2019-2024', ylab='Maternal Deaths',xlab='Time',col=colors)
legend("topleft",                     
       legend = c("Hispanic", "Black, Non-Hispanic", "White, Non-Hispanic",'Asian, Non-Hispanic','American Indian, Alaska Native, Non-Hisp'),  # Labels
       col = colors,                   
       lty = 1,                        
       lwd = 2,                     
       bty = "n",
       cex = 0.8)                      

HisR <- filtered_df%>%
  filter(Subgroup=='Hispanic')%>%
  select(Maternal_Mortality_Rate)%>%
  rename(Hispanic=Maternal_Mortality_Rate)
blR <- filtered_df%>%
  filter(Subgroup=='Black, Non-Hispanic')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('Black_Non-Hispanic' = Maternal_Mortality_Rate)
whR <- filtered_df%>%
  filter(Subgroup=='White, Non-Hispanic')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('White_Non-Hispanic'=Maternal_Mortality_Rate)
asR <- filtered_df%>%
  filter(Subgroup=='Asian, Non-Hispanic')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('Asian_Non-Hispanic'=Maternal_Mortality_Rate)
ainaR <- df%>%
  filter(Subgroup=='American Indian or Alaska Native, Non-Hispanic')%>%
  select(Maternal_Mortality_Rate)%>% 
  rename('American_Indian_or_Alaska_Native, Non-Hispanic'=Maternal_Mortality_Rate)
ethnR_ts <- bind_cols(HisR,blR,whR,asR,ainaR) 
ethnR_ts <- ts(ethnR_ts, start=c(2019,1), frequency = 12)
colors=c('red','blue','green','orange','yellow')
plot.ts(ethnR_ts, plot.type='single',main='Monthly Maternal Mortality Rate by Ethnicity 2019-2024', ylab='Maternal Deaths per 100,000 Live Births',xlab='Time',col=colors)
legend("topleft",                     
       legend = c("Hispanic", "Black, Non-Hispanic", "White, Non-Hispanic",'Asian, Non-Hispanic','American Indian, Alaska Native, Non-Hisp'),  # Labels
       col = colors,                   
       lty = 1,                        
       lwd = 2,                     
       bty = "n",
       cex = 0.7)            

ageu25 <- deaths_df%>%
  filter(Subgroup=='Under 25 years')%>%
  select(Maternal_Deaths)%>%
  rename('Under_25_years'=Maternal_Deaths)
age25 <- filtered_df%>%
  filter(Subgroup=='25-39 years')%>%
  select(Maternal_Deaths)%>%
  rename('25-39_years' = Maternal_Deaths)
ageo40 <- filtered_df%>%
  filter(Subgroup=='40 years and over')%>%
  select(Maternal_Deaths)%>%
  rename('40_years_and_over'=Maternal_Deaths)
age_ts <- bind_cols(ageu25,age25,ageo40) 
age_ts <- ts(age_ts, start=c(2019,1), frequency = 12)
colors=c('red','blue','green')
plot.ts(age_ts, plot.type='single',main='Monthly Maternal Deaths by Age Group 2019-2024', ylab='Maternal Deaths',xlab='Time',col=colors)
legend("topleft",                     
       legend = c("Under 25 Years", "25-39 Years", "40 Years and over"),  # Labels
       col = colors,                   
       lty = 1,                        
       lwd = 2,                     
       bty = "n",
       cex = 0.8) 

ageu25R <- filtered_df%>%
  filter(Subgroup=='Under 25 years')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('Under_25_years'=Maternal_Mortality_Rate)
age25R <- filtered_df%>%
  filter(Subgroup=='25-39 years')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('25-39_years' = Maternal_Mortality_Rate)
ageo40R <- filtered_df%>%
  filter(Subgroup=='40 years and over')%>%
  select(Maternal_Mortality_Rate)%>%
  rename('40_years_and_over'=Maternal_Mortality_Rate)
ageR_ts <- bind_cols(ageu25R,age25R,ageo40R) 
ageR_ts <- ts(ageR_ts, start=c(2019,1), frequency = 12)
colors=c('red','blue','green')
plot.ts(ageR_ts, plot.type='single',main='Monthly Maternal Mortality Rate by Age Group 2019-2024', ylab='Maternal Deaths per 100,000 live Births',xlab='Time',col=colors)
legend("topleft",                     
       legend = c("Under 25 Years", "25-39 Years", "40 Years and over"),  # Labels
       col = colors,                   
       lty = 1,                        
       lwd = 2,                     
       bty = "n",
       cex = 0.8) 

plot_histogram(deaths_df_ethnic$Maternal_Deaths,title='Histogram of Maternal Deaths - Ethnicity')

plot_histogram(deaths_df_age$Maternal_Deaths, title='Histogram of Maternal Deaths - Age')

plot_histogram(total_df$Maternal_Deaths,title='Histogram of Maternal Deaths - Total')

#check overdispersion by for ethnicity data
deaths_df_ethnic %>% summarize(mean(Maternal_Deaths),
var(Maternal_Deaths))

  mean(Maternal_Deaths) var(Maternal_Deaths)
1              165.5779              18843.1

#check overdispersion for age data
deaths_df_age %>% summarize(mean(Maternal_Deaths),var(Maternal_Deaths))

  mean(Maternal_Deaths) var(Maternal_Deaths)
1              276.2407              54102.1

#check overdispersion for total data
total_df %>% summarize(mean(Maternal_Deaths),
var(Maternal_Deaths))

  mean(Maternal_Deaths) var(Maternal_Deaths)
1              828.8889             30819.62

<<<<<<< HEAD

5ff3583d31d40c93aaa9005a06c34bf67126bfd1

poissonmodel_glmm_all <- glmer(Maternal_Deaths ~ Ethnicity +Age_Group+ Dobbs_Era+(1|Year),
                              offset=log(Live_Births),
                              family = poisson(link = "log"), 
                              data = deaths_df3)
poissonmodel_glmm_all

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: Maternal_Deaths ~ Ethnicity + Age_Group + Dobbs_Era + (1 | Year)
   Data: deaths_df3
 Offset: log(Live_Births)
      AIC       BIC    logLik -2*log(L)  df.resid 
 5133.339  5176.778 -2556.670  5113.339       559 
Random effects:
 Groups Name        Std.Dev.
 Year   (Intercept) 0.1755  
Number of obs: 569, groups:  Year, 6
Fixed Effects:
                                            (Intercept)  
                                                -8.7966  
                           EthnicityBlack, Non-Hispanic  
                                                 1.3112  
                           EthnicityWhite, Non-Hispanic  
                                                 0.2857  
                                      EthnicityHispanic  
                                                 0.1735  
EthnicityAmerican Indian or Alaska Native, Non-Hispanic  
                                                 1.4393  
                                       EthnicityUnknown  
                                                 0.0412  
                                   Age_Group25-39 years  
                                                 0.3926  
                             Age_Group40 years and over  
                                                 1.8197  
                                    Dobbs_EraPost-Dobbs  
                                                -0.2274  
fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

check_overdispersion(poissonmodel_glmm_all)

# Overdispersion test

       dispersion ratio =    2.318
  Pearson's Chi-Squared = 1295.809
                p-value =  < 0.001

##overdispersion still detected

install.packages("gglmmTMB")

##with offset # all fixed effects
all_glmmodel_nb <- glmmTMB(
  Maternal_Deaths ~  Ethnicity + Age_Group +Dobbs_Era+ (1|Year),
  offset=log(Live_Births),
  family = nbinom2,
  data = deaths_df3
)
summary(all_glmmodel_nb)

 Family: nbinom2  ( log )
Formula:          
Maternal_Deaths ~ Ethnicity + Age_Group + Dobbs_Era + (1 | Year)
Data: deaths_df3
 Offset: log(Live_Births)

      AIC       BIC    logLik -2*log(L)  df.resid 
   4822.7    4870.5   -2400.4    4800.7       558 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Year   (Intercept) 0.03253  0.1804  
Number of obs: 569, groups:  Year, 6

Dispersion parameter for nbinom2 family ():  149 

Conditional model:
                                                        Estimate Std. Error
(Intercept)                                             -8.79484    0.07779
EthnicityBlack, Non-Hispanic                             1.31589    0.02624
EthnicityWhite, Non-Hispanic                             0.28350    0.02601
EthnicityHispanic                                        0.16211    0.02701
EthnicityAmerican Indian or Alaska Native, Non-Hispanic  1.43374    0.04033
EthnicityUnknown                                         0.04376    0.02747
Age_Group25-39 years                                     0.38833    0.01816
Age_Group40 years and over                               1.81026    0.02049
Age_GroupUnknown                                              NA         NA
Dobbs_EraPost-Dobbs                                     -0.22654    0.02287
                                                        z value Pr(>|z|)    
(Intercept)                                             -113.06  < 2e-16 ***
EthnicityBlack, Non-Hispanic                              50.16  < 2e-16 ***
EthnicityWhite, Non-Hispanic                              10.90  < 2e-16 ***
EthnicityHispanic                                          6.00 1.94e-09 ***
EthnicityAmerican Indian or Alaska Native, Non-Hispanic   35.55  < 2e-16 ***
EthnicityUnknown                                           1.59    0.111    
Age_Group25-39 years                                      21.38  < 2e-16 ***
Age_Group40 years and over                                88.35  < 2e-16 ***
Age_GroupUnknown                                             NA       NA    
Dobbs_EraPost-Dobbs                                       -9.91  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ethnicity_agegroup_glmmodel_nb <- glmmTMB(
  Maternal_Deaths ~  Ethnicity+ Age_Group + (1|Year),
  offset=log(Live_Births),
  family = nbinom2,
  data = deaths_df3
)
summary(ethnicity_agegroup_glmmodel_nb)

 Family: nbinom2  ( log )
Formula:          Maternal_Deaths ~ Ethnicity + Age_Group + (1 | Year)
Data: deaths_df3
 Offset: log(Live_Births)

      AIC       BIC    logLik -2*log(L)  df.resid 
   4908.8    4952.2   -2444.4    4888.8       559 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Year   (Intercept) 0.03363  0.1834  
Number of obs: 569, groups:  Year, 6

Dispersion parameter for nbinom2 family ():  107 

Conditional model:
                                                        Estimate Std. Error
(Intercept)                                             -8.88758    0.07864
EthnicityBlack, Non-Hispanic                             1.31609    0.02760
EthnicityWhite, Non-Hispanic                             0.28274    0.02739
EthnicityHispanic                                        0.16064    0.02833
EthnicityAmerican Indian or Alaska Native, Non-Hispanic  1.43282    0.04129
EthnicityUnknown                                         0.04379    0.02878
Age_Group25-39 years                                     0.38801    0.02008
Age_Group40 years and over                               1.80885    0.02221
Age_GroupUnknown                                              NA         NA
                                                        z value Pr(>|z|)    
(Intercept)                                             -113.02  < 2e-16 ***
EthnicityBlack, Non-Hispanic                              47.68  < 2e-16 ***
EthnicityWhite, Non-Hispanic                              10.32  < 2e-16 ***
EthnicityHispanic                                          5.67 1.43e-08 ***
EthnicityAmerican Indian or Alaska Native, Non-Hispanic   34.70  < 2e-16 ***
EthnicityUnknown                                           1.52    0.128    
Age_Group25-39 years                                      19.33  < 2e-16 ***
Age_Group40 years and over                                81.45  < 2e-16 ***
Age_GroupUnknown                                             NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##with no offset # all fixed effects
allno_glmmodel_nb <- glmmTMB(
  Maternal_Deaths ~  Ethnicity + Age_Group+ Dobbs_Era + (1 |Year),
  family = nbinom2,
  data = deaths_df3
)
summary(allno_glmmodel_nb)

 Family: nbinom2  ( log )
Formula:          
Maternal_Deaths ~ Ethnicity + Age_Group + Dobbs_Era + (1 | Year)
Data: deaths_df3

      AIC       BIC    logLik -2*log(L)  df.resid 
   4818.1    4865.9   -2398.1    4796.1       558 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Year   (Intercept) 0.02963  0.1721  
Number of obs: 569, groups:  Year, 6

Dispersion parameter for nbinom2 family ():  156 

Conditional model:
                                                        Estimate Std. Error
(Intercept)                                              3.51807    0.07458
EthnicityBlack, Non-Hispanic                             2.16074    0.02608
EthnicityWhite, Non-Hispanic                             2.40614    0.02586
EthnicityHispanic                                        1.56987    0.02686
EthnicityAmerican Indian or Alaska Native, Non-Hispanic -0.70147    0.04022
EthnicityUnknown                                         1.33578    0.02732
Age_Group25-39 years                                     1.59882    0.01793
Age_Group40 years and over                               0.03840    0.02028
Age_GroupUnknown                                              NA         NA
Dobbs_EraPost-Dobbs                                     -0.22821    0.02250
                                                        z value Pr(>|z|)    
(Intercept)                                               47.17   <2e-16 ***
EthnicityBlack, Non-Hispanic                              82.85   <2e-16 ***
EthnicityWhite, Non-Hispanic                              93.06   <2e-16 ***
EthnicityHispanic                                         58.45   <2e-16 ***
EthnicityAmerican Indian or Alaska Native, Non-Hispanic  -17.44   <2e-16 ***
EthnicityUnknown                                          48.89   <2e-16 ***
Age_Group25-39 years                                      89.18   <2e-16 ***
Age_Group40 years and over                                 1.89   0.0583 .  
Age_GroupUnknown                                             NA       NA    
Dobbs_EraPost-Dobbs                                      -10.14   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ethnicity_agegroupno_glmmodel_nb <- glmmTMB(
  Maternal_Deaths ~  Ethnicity+ Age_Group + (1 |Year),
  family = nbinom2,
  data = deaths_df3
)
summary(ethnicity_agegroupno_glmmodel_nb)

 Family: nbinom2  ( log )
Formula:          Maternal_Deaths ~ Ethnicity + Age_Group + (1 | Year)
Data: deaths_df3

      AIC       BIC    logLik -2*log(L)  df.resid 
   4907.7    4951.1   -2443.9    4887.7       559 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Year   (Intercept) 0.03167  0.178   
Number of obs: 569, groups:  Year, 6

Dispersion parameter for nbinom2 family ():  110 

Conditional model:
                                                        Estimate Std. Error
(Intercept)                                              3.42510    0.07652
EthnicityBlack, Non-Hispanic                             2.16034    0.02749
EthnicityWhite, Non-Hispanic                             2.40499    0.02727
EthnicityHispanic                                        1.56804    0.02823
EthnicityAmerican Indian or Alaska Native, Non-Hispanic -0.70298    0.04121
EthnicityUnknown                                         1.33526    0.02867
Age_Group25-39 years                                     1.59860    0.01991
Age_Group40 years and over                               0.03684    0.02205
Age_GroupUnknown                                              NA         NA
                                                        z value Pr(>|z|)    
(Intercept)                                               44.76   <2e-16 ***
EthnicityBlack, Non-Hispanic                              78.59   <2e-16 ***
EthnicityWhite, Non-Hispanic                              88.18   <2e-16 ***
EthnicityHispanic                                         55.55   <2e-16 ***
EthnicityAmerican Indian or Alaska Native, Non-Hispanic  -17.06   <2e-16 ***
EthnicityUnknown                                          46.57   <2e-16 ***
Age_Group25-39 years                                      80.30   <2e-16 ***
Age_Group40 years and over                                 1.67   0.0948 .  
Age_GroupUnknown                                             NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

performance::check_model(all_glmmodel_nb, verbose = TRUE,residual_type = "normal")#interesting

performance::r2(all_glmmodel_nb)

# R2 for Mixed Models

  Conditional R2: 0.060
     Marginal R2: 0.056

check_singularity(all_glmmodel_nb)#random effect is not a singularity

[1] FALSE

check_collinearity(all_glmmodel_nb)

# Check for Multicollinearity

Low Correlation

      Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 Ethnicity 2.57 [2.27,     2.94]     1.60      0.39     [0.34, 0.44]
 Age_Group 2.57 [2.27,     2.94]     1.60      0.39     [0.34, 0.44]
 Dobbs_Era 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

performance::check_model(ethnicity_agegroup_glmmodel_nb, verbose=TRUE, residual_type = "normal")#interesting

performance::r2(ethnicity_agegroup_glmmodel_nb)

# R2 for Mixed Models

  Conditional R2: 0.059
     Marginal R2: 0.055

check_singularity(ethnicity_agegroup_glmmodel_nb)

[1] FALSE

check_collinearity(ethnicity_agegroup_glmmodel_nb)

# Check for Multicollinearity

Low Correlation

      Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 Ethnicity 2.47 [2.19, 2.83]     1.57      0.40     [0.35, 0.46]
 Age_Group 2.47 [2.19, 2.83]     1.57      0.40     [0.35, 0.46]

performance::check_model(allno_glmmodel_nb, verbose=TRUE, residual_type = "normal")##interesting

performance::r2(allno_glmmodel_nb)

# R2 for Mixed Models

  Conditional R2: 0.991
     Marginal R2: 0.969

check_singularity(allno_glmmodel_nb)

[1] FALSE

check_collinearity(allno_glmmodel_nb)

# Check for Multicollinearity

Low Correlation

      Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 Ethnicity 2.59 [2.29,     2.96]     1.61      0.39     [0.34, 0.44]
 Age_Group 2.59 [2.29,     2.96]     1.61      0.39     [0.34, 0.44]
 Dobbs_Era 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

performance::check_model(ethnicity_agegroupno_glmmodel_nb, verbose=TRUE, residual_type = "normal")#interesting

performance::r2(ethnicity_agegroupno_glmmodel_nb)

# R2 for Mixed Models

  Conditional R2: 0.989
     Marginal R2: 0.965

check_singularity(ethnicity_agegroupno_glmmodel_nb)

[1] FALSE

check_collinearity(ethnicity_agegroupno_glmmodel_nb)

# Check for Multicollinearity

Low Correlation

      Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 Ethnicity 2.48 [2.19, 2.84]     1.58      0.40     [0.35, 0.46]
 Age_Group 2.48 [2.19, 2.84]     1.58      0.40     [0.35, 0.46]

resid_pearson <- residuals(all_glmmodel_nb, type = "pearson")
plot(resid_pearson ~ fitted(all_glmmodel_nb), main = "Pearson residuals vs Fitted")
abline(h = 0, col = "red")

overdisp_fun <- function(all_glmmodel_nb) {
  rdf <- df.residual(all_glmmodel_nb)
  rp <- residuals(all_glmmodel_nb, type = "pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq / rdf
  pval <- pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE)
  c(chisq = Pearson.chisq, ratio = prat, rdf = rdf, p = pval)
}

overdisp_fun(all_glmmodel_nb)

      chisq       ratio         rdf           p 
553.4165496   0.9917859 558.0000000   0.5468108 

#ratio is around 1 which is good and not too indicative of overdispersion
#pp-value is greater that P<0.005 so that means there is no significant overdispersion

arm::binnedplot(fitted(all_glmmodel_nb), 
           residuals(all_glmmodel_nb, type = "response"), 
           main = "Binned residual plot")

# fan shape  heteroskedasticity (non-constant variance)- coon when with count data

ranef(all_glmmodel_nb)$cond

$Year
      (Intercept)
2019 -0.263523739
2020 -0.142158611
2021  0.125467718
2022  0.297797811
2023 -0.013392036
2024 -0.004252477

re <- ranef(all_glmmodel_nb)
cond_re <- re$cond
year_re <- as.data.frame(cond_re$Year)
year_re$Year <- rownames(year_re)

ggplot(year_re, aes(x = factor(Year), y = `(Intercept)`)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Random Effect Intercepts by Year",
       x = "Year",
       y = "Random Effect Estimate (Intercept)") +
  theme_minimal()

#dashed zero is the overall random intercept
# points above zero - years with above average random effect
#points below zero- years below average random effect

sjPlot::plot_model(all_glmmodel_nb, type = "re") #no extreme points with random effects

sjPlot::tab_model(all_glmmodel_nb)

	Maternal Deaths
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	0.00	0.00 – 0.00	<0.001
Ethnicity [Black, Non-Hispanic]	3.73	3.54 – 3.92	<0.001
Ethnicity [White, Non-Hispanic]	1.33	1.26 – 1.40	<0.001
Ethnicity [Hispanic]	1.18	1.12 – 1.24	<0.001
Ethnicity [American Indian or Alaska Native, Non-Hispanic]	4.19	3.88 – 4.54	<0.001
Ethnicity [Unknown]	1.04	0.99 – 1.10	0.111
Age_Group25-39 years	1.47	1.42 – 1.53	<0.001
Age Group [40 years and over]	6.11	5.87 – 6.36	<0.001
Dobbs Era [Post-Dobbs]	0.80	0.76 – 0.83	<0.001
Random Effects
σ2	7.96
τ00 Year	0.03
ICC	0.00
N Year	6
Observations	569
Marginal R2 / Conditional R2	0.056 / 0.060

resid_pearson <- residuals(ethnicity_agegroup_glmmodel_nb, type = "pearson")
plot(resid_pearson ~ fitted(ethnicity_agegroup_glmmodel_nb), main = "Pearson residuals vs Fitted")
abline(h = 0, col = "red")

overdisp_fun <- function(ethnicity_agegroup_glmmodel_nb) {
  rdf <- df.residual(ethnicity_agegroup_glmmodel_nb)
  rp <- residuals(ethnicity_agegroup_glmmodel_nb, type = "pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq / rdf
  pval <- pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE)
  c(chisq = Pearson.chisq, ratio = prat, rdf = rdf, p = pval)
}

overdisp_fun(ethnicity_agegroup_glmmodel_nb)

      chisq       ratio         rdf           p 
538.2717872   0.9629191 559.0000000   0.7283260 

#ratio is around 1 which is good and not too indicative of overdispersion
#pp-value is greater that P<0.005 so that means there is no significant overdispersion

arm::binnedplot(fitted(ethnicity_agegroup_glmmodel_nb), 
           residuals(ethnicity_agegroup_glmmodel_nb, type = "response"), 
           main = "Binned residual plot")

# fan shape  heteroskedasticity (non-constant variance)- coon when with count data

ranef(ethnicity_agegroup_glmmodel_nb)$cond

$Year
     (Intercept)
2019 -0.17024416
2020 -0.05111811
2021  0.21881739
2022  0.28491518
2023 -0.14883772
2024 -0.13365740

re <- ranef(ethnicity_agegroup_glmmodel_nb)
cond_re <- re$cond
year_re <- as.data.frame(cond_re$Year)
year_re$Year <- rownames(year_re)

ggplot(year_re, aes(x = factor(Year), y = `(Intercept)`)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Random Effect Intercepts by Year",
       x = "Year",
       y = "Random Effect Estimate (Intercept)") +
  theme_minimal()

#dashed zero is the overall random intercept
# points above zero - years with above average random effect
#points below zero- years below average random effect

sjPlot::plot_model(ethnicity_agegroup_glmmodel_nb, type = "re") #no extreme points with random effects

sjPlot::tab_model(ethnicity_agegroup_glmmodel_nb)

	Maternal Deaths
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	0.00	0.00 – 0.00	<0.001
Ethnicity [Black, Non-Hispanic]	3.73	3.53 – 3.94	<0.001
Ethnicity [White, Non-Hispanic]	1.33	1.26 – 1.40	<0.001
Ethnicity [Hispanic]	1.17	1.11 – 1.24	<0.001
Ethnicity [American Indian or Alaska Native, Non-Hispanic]	4.19	3.86 – 4.54	<0.001
Ethnicity [Unknown]	1.04	0.99 – 1.11	0.128
Age_Group25-39 years	1.47	1.42 – 1.53	<0.001
Age Group [40 years and over]	6.10	5.84 – 6.37	<0.001
Random Effects
σ2	7.96
τ00 Year	0.03
ICC	0.00
N Year	6
Observations	569
Marginal R2 / Conditional R2	0.055 / 0.059

resid_pearson <- residuals(allno_glmmodel_nb, type = "pearson")
plot(resid_pearson ~ fitted(allno_glmmodel_nb), main = "Pearson residuals vs Fitted")
abline(h = 0, col = "red")

overdisp_fun <- function(allno_glmmodel_nb) {
  rdf <- df.residual(allno_glmmodel_nb)
  rp <- residuals(allno_glmmodel_nb, type = "pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq / rdf
  pval <- pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE)
  c(chisq = Pearson.chisq, ratio = prat, rdf = rdf, p = pval)
}

overdisp_fun(allno_glmmodel_nb)

      chisq       ratio         rdf           p 
558.3912206   1.0007011 558.0000000   0.4873697 

#ratio is around 1 which is good and not too indicative of overdispersion
#pp-value is greater that P<0.005 so that means there is no significant overdispersion

arm::binnedplot(fitted(allno_glmmodel_nb), 
           residuals(allno_glmmodel_nb, type = "response"), 
           main = "Binned residual plot")

# fan shape  heteroskedasticity (non-constant variance)- coon when with count data

ranef(allno_glmmodel_nb)$cond

$Year
     (Intercept)
2019 -0.23906046
2020 -0.13342844
2021  0.10309415
2022  0.30134973
2023 -0.01548869
2024 -0.01651757

re <- ranef(allno_glmmodel_nb)
cond_re <- re$cond
year_re <- as.data.frame(cond_re$Year)
year_re$Year <- rownames(year_re)

ggplot(year_re, aes(x = factor(Year), y = `(Intercept)`)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Random Effect Intercepts by Year",
       x = "Year",
       y = "Random Effect Estimate (Intercept)") +
  theme_minimal()

#dashed zero is the overall random intercept
# points above zero - years with above average random effect
#points below zero- years below average random effect

sjPlot::tab_model(allno_glmmodel_nb)#no offset

	Maternal Deaths
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	33.72	29.13 – 39.03	<0.001
Ethnicity [Black, Non-Hispanic]	8.68	8.25 – 9.13	<0.001
Ethnicity [White, Non-Hispanic]	11.09	10.54 – 11.67	<0.001
Ethnicity [Hispanic]	4.81	4.56 – 5.07	<0.001
Ethnicity [American Indian or Alaska Native, Non-Hispanic]	0.50	0.46 – 0.54	<0.001
Ethnicity [Unknown]	3.80	3.60 – 4.01	<0.001
Age_Group25-39 years	4.95	4.78 – 5.12	<0.001
Age Group [40 years and over]	1.04	1.00 – 1.08	0.058
Dobbs Era [Post-Dobbs]	0.80	0.76 – 0.83	<0.001
Random Effects
σ2	0.01
τ00 Year	0.03
ICC	0.73
N Year	6
Observations	569
Marginal R2 / Conditional R2	0.969 / 0.991

resid_pearson <- residuals(ethnicity_agegroupno_glmmodel_nb, type = "pearson")
plot(resid_pearson ~ fitted(ethnicity_agegroupno_glmmodel_nb), main = "Pearson residuals vs Fitted")
abline(h = 0, col = "red")

overdisp_fun <- function(ethnicity_agegroupno_glmmodel_nb) {
  rdf <- df.residual(ethnicity_agegroupno_glmmodel_nb)
  rp <- residuals(ethnicity_agegroupno_glmmodel_nb, type = "pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq / rdf
  pval <- pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE)
  c(chisq = Pearson.chisq, ratio = prat, rdf = rdf, p = pval)
}

overdisp_fun(ethnicity_agegroupno_glmmodel_nb)

      chisq       ratio         rdf           p 
542.2582225   0.9700505 559.0000000   0.6864459 

#ratio is around 1 which is good and not too indicative of overdispersion
#pp-value is greater that P<0.005 so that means there is no significant overdispersion

arm::binnedplot(fitted(ethnicity_agegroupno_glmmodel_nb), 
           residuals(ethnicity_agegroupno_glmmodel_nb, type = "response"), 
           main = "Binned residual plot")

# fan shape  heteroskedasticity (non-constant variance)- coon when with count data

ranef(ethnicity_agegroupno_glmmodel_nb)$cond

$Year
     (Intercept)
2019 -0.14482143
2020 -0.04174171
2021  0.19656817
2022  0.28820781
2023 -0.15191772
2024 -0.14644654

re <- ranef(ethnicity_agegroupno_glmmodel_nb)
cond_re <- re$cond
year_re <- as.data.frame(cond_re$Year)
year_re$Year <- rownames(year_re)

ggplot(year_re, aes(x = factor(Year), y = `(Intercept)`)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Random Effect Intercepts by Year",
       x = "Year",
       y = "Random Effect Estimate (Intercept)") +
  theme_minimal()

#dashed zero is the overall random intercept
# points above zero - years with above average random effect
#points below zero- years below average random effect

sjPlot::tab_model(ethnicity_agegroupno_glmmodel_nb)

	Maternal Deaths
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	30.73	26.45 – 35.70	<0.001
Ethnicity [Black, Non-Hispanic]	8.67	8.22 – 9.15	<0.001
Ethnicity [White, Non-Hispanic]	11.08	10.50 – 11.69	<0.001
Ethnicity [Hispanic]	4.80	4.54 – 5.07	<0.001
Ethnicity [American Indian or Alaska Native, Non-Hispanic]	0.50	0.46 – 0.54	<0.001
Ethnicity [Unknown]	3.80	3.59 – 4.02	<0.001
Age_Group25-39 years	4.95	4.76 – 5.14	<0.001
Age Group [40 years and over]	1.04	0.99 – 1.08	0.095
Random Effects
σ2	0.01
τ00 Year	0.03
ICC	0.70
N Year	6
Observations	569
Marginal R2 / Conditional R2	0.965 / 0.989

##Comare AIC and see who has influential points and then we find best model

Modeling and Results

• Explain your data preprocessing and cleaning steps.

• Present your key findings in a clear and concise manner.

• Use visuals to support your claims.

• Tell a story about what the data reveals.

Conclusion

• Summarize your key findings.

• Discuss the implications of your results.

References

Agresti, A. 2015. Foundations of Linear and Generalized Linear Models. Wiley Series in Probability and Statistics. Wiley. https://books.google.com/books?id=jlIqBgAAQBAJ.

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Candy, Steven G. 2000. “The Application of Generalized Linear Mixed Models to Multi-Level Sampling for Insect Population Monitoring.” Environmental and Ecological Statistics 7 (3): 217–38.

Lee, Youngjo, and John A Nelder. 1996. “Hierarchical Generalized Linear Models.” Journal of the Royal Statistical Society Series B: Statistical Methodology 58 (4): 619–56.

McGillycuddy, Maeve, David I. Warton, Gordana Popovic, and Benjamin M. Bolker. 2025. “Parsimoniously Fitting Large Multivariate Random Effects in glmmTMB.” Journal of Statistical Software 112 (1): 1–19. https://doi.org/10.18637/jss.v112.i01.

Owili, Patrick Opiyo, Tang-Huang Lin, Miriam Adoyo Muga, and Wei-Hung Lien. 2020. “Impacts of Discriminated PM2. 5 on Global Under-Five and Maternal Mortality.” Scientific Reports 10 (1): 17654.

R Core Team. 2025. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Salinas Ruı́z, Josafhat, Osval Antonio Montesinos López, Gabriela Hernández Ramı́rez, and Jose Crossa Hiriart. 2023. Generalized Linear Mixed Models with Applications in Agriculture and Biology. Springer Nature.

Tawiah, Kassim, Samuel Iddi, and Anani Lotsi. 2020. “On Zero-Inflated Hierarchical Poisson Models with Application to Maternal Mortality Data.” International Journal of Mathematics and Mathematical Sciences 2020 (1): 1407320.

Thall, Peter F. 1988. “Mixed Poisson Likelihood Regression Models for Longitudinal Interval Count Data.” Biometrics, 197–209.

Thompson, Jennifer A, Clemence Leyrat, Katherine L Fielding, and Richard J Hayes. 2022. “Cluster Randomised Trials with a Binary Outcome and a Small Number of Clusters: Comparison of Individual and Cluster Level Analysis Method.” BMC Medical Research Methodology 22 (1): 222.

Wang, Ke-Sheng, Xuefeng Liu, Muyiwa Ategbole, Xin Xie, Ying Liu, Chun Xu, Changchun Xie, and Zhanxin Sha. 2017. “Generalized Linear Mixed Model Analysis of Urban-Rural Differences in Social and Behavioral Factors for Colorectal Cancer Screening.” Asian Pacific Journal of Cancer Prevention: APJCP 18 (9): 2581.

Zuur, Alain F, Elena N Ieno, Neil J Walker, Anatoly A Saveliev, Graham M Smith, et al. 2009. Mixed Effects Models and Extensions in Ecology with r. Vol. 574. Springer.